#setwd('Dropbox/Mycobiome Xav Analysis/')
######### LIBRARIES
library(phyloseq)
library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.0 ──
## ✓ ggplot2 3.3.3 ✓ purrr 0.3.4
## ✓ tibble 3.1.0 ✓ dplyr 1.0.5
## ✓ tidyr 1.1.3 ✓ stringr 1.4.0
## ✓ readr 1.4.0 ✓ forcats 0.5.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
library(vegan)
## Loading required package: permute
## Loading required package: lattice
## This is vegan 2.5-7
library(brms)
## Loading required package: Rcpp
## Loading 'brms' package (version 2.15.0). Useful instructions
## can be found by typing help('brms'). A more detailed introduction
## to the package is available through vignette('brms_overview').
##
## Attaching package: 'brms'
## The following object is masked from 'package:phyloseq':
##
## nsamples
## The following object is masked from 'package:stats':
##
## ar
library(dplyr)
library(tidybayes)
##
## Attaching package: 'tidybayes'
## The following objects are masked from 'package:brms':
##
## dstudent_t, pstudent_t, qstudent_t, rstudent_t
library(tidyr)
library(cowplot)
library(microbiome)
##
## microbiome R package (microbiome.github.com)
##
##
##
## Copyright (C) 2011-2020 Leo Lahti,
## Sudarshan Shetty et al. <microbiome.github.io>
##
## Attaching package: 'microbiome'
## The following object is masked from 'package:vegan':
##
## diversity
## The following object is masked from 'package:ggplot2':
##
## alpha
## The following object is masked from 'package:base':
##
## transform
options(mc.cores=parallel::detectCores())
#Load Data Files
myco<-readRDS('clean.its.rds')
micro<-readRDS('clean.16s.rds')
##Metadata
metadata<-read.csv('Sample Metadata with Diet May20.csv',header=T)
#head(metadata)
##RNG SEED
set.seed(21042020)
## Use multiple Cores for Model Fitting (brms)
options(mc.cores=parallel::detectCores())
#Global Plot Options
plotopts<- theme(axis.text=element_text(size=12),axis.title=element_text(size=16),strip.text=element_text(size=20),legend.title = element_text(size=18),legend.text = element_text(size=18))
#library(devtools)
#install_github("zdk123/SpiecEasi")
### Core Datasets of Myco and Micro Sample
microbial_common_samples<-Reduce(intersect, list(rownames(sample_data(myco)),rownames(sample_data(micro))))
myco_common<-subset_samples(myco,Sample %in% microbial_common_samples)
micro_common<-subset_samples(micro,Sample %in% microbial_common_samples)
#Checks
phyloseq::nsamples(myco_common) == phyloseq::nsamples(micro_common)
## [1] TRUE
all.equal(as.character(sample_data(myco_common)$Sample),as.character(sample_data(micro_common)$Sample)) #TRUE
## [1] TRUE
####### MERGING OBJECTS
#rename if required and check again
taxa_names(micro_common) <- paste0("bacteria-", taxa_names(micro_common))
#taxa_names(micro_common)
taxa_names(myco_common) <- paste0("fungi-", taxa_names(myco_common))
#taxa_names(myco_common)
#Make Sure No Duplicate SV Names
table(taxa_names(micro_common) %in% taxa_names(myco_common)) ## ALL FALSE
##
## FALSE
## 803
#Merge Objects
microbial_merge<-merge_phyloseq(micro_common,myco_common)
#More Checking - Samples and Microbial Taxa Totals
ntaxa(microbial_merge) == (ntaxa(micro_common) + ntaxa(myco_common)) #TRUE
## [1] TRUE
phyloseq::nsamples(microbial_merge) == length(microbial_common_samples) #TRUE
## [1] TRUE
### Metadata Checking
table(sample_data(microbial_merge)$Common_name)
##
## African palm weevil larvae Bank vole
## 6 10
## Barred hamlet Blackcap
## 5 5
## Blue tit Brown shrimp
## 3 2
## Capuchin monkey Carrion crow
## 9 2
## Chiffchaff Cockroach
## 3 7
## Collared dove Common midwife toad
## 5 11
## Common toad Eastern black rhino
## 2 10
## European eel Foureye butterflyfish
## 1 9
## Goldfinch Great tit
## 3 3
## Greater white-toothed shrew Grey squirrel
## 7 5
## Hard tick Hedgehog
## 2 11
## Honey bee Lesser horseshoe bat
## 6 2
## Light-bellied brent goose Northern muriqui
## 9 10
## Phofung river frog Pygmy shrew
## 4 10
## Red deer Red squirrel
## 10 11
## Reed bunting Reed warbler
## 2 2
## Robin Roe deer
## 4 4
## Song thrush Stock dove
## 5 1
## Striped field mouse Tsetse fly
## 10 6
## Turtle dove Vase sponge
## 7 10
## Wild pony Wood mouse
## 10 9
## Woodpigeon Yellow-necked mouse
## 4 10
## Yellowhammer Yellowhead wrasse
## 5 5
########## SPIEC-EASI MODELSs
library(SpiecEasi)
library(igraph)
### Rename Phyla
tax_table(microbial_merge)[,2][which(tax_table(microbial_merge)[,2]=="Kingdom_k__Fungi")] <-"Kingdom_Fungi"
tax_table(microbial_merge)[,2]<-gsub("^p__","",tax_table(microbial_merge)[,2])
## Function for Colour Edges by Sign of Relationship
extract_edge_cols<-function(physeq_obj,speic_obj,igraph_obj){
otu.ids=rownames(tax_table(physeq_obj))
betaMat=as.matrix(symBeta(getOptBeta(speic_obj)))
edges=E(igraph_obj)
edge.colors=c()
for(e.index in 1:length(edges)){
adj.nodes=ends(igraph_obj,edges[e.index])
xindex=which(otu.ids==adj.nodes[1])
yindex=which(otu.ids==adj.nodes[2])
beta=betaMat[xindex,yindex]
if(beta>0){
edge.colors=append(edge.colors,"forestgreen")
}else if(beta<0){
edge.colors=append(edge.colors,"red")
}
}
return (edge.colors)
}
### Palette for Kingdom
library(RColorBrewer)
length(unique(tax_table(microbial_merge)[,2])) #21 Phyla
## [1] 21
global_phy_cols<-data.frame(Phylum=unique(tax_table(microbial_merge)[,2]))
#global_phy_cols$color<-colorRampPalette(brewer.pal(11,"Spectral"))(nrow(global_phy_cols))
global_phy_cols$color<-c(brewer.pal(12,"Paired"),brewer.pal(9,"Pastel1"))
### Mammal Test
mammals <- subset_samples(microbial_merge,Taxa=="Mammal") %>% filter_taxa(function(x) sum(x > 3) > (0.05*length(x)), TRUE)
mammal_speic<- mammals %>% spiec.easi(method='mb', lambda.min.ratio=1e-2, nlambda=20,pulsar.params=list(rep.num=50))
## Applying data transformations...
## Selecting model with pulsar using stars...
## Fitting final estimate with mb...
## done
ig.mammal <- adj2igraph(getRefit(mammal_speic),vertex.attr=list(name=taxa_names(mammals),Kingdom=tax_table(mammals)[,1],Phylum=tax_table(mammals)[,2]))
#Assign Shapes to Kingdoms + Colours to Phylum
V(ig.mammal)$shape<-ifelse(V(ig.mammal)$Kingdom=="Bacteria","circle","square")
V(ig.mammal)$color<-global_phy_cols$color[match(V(ig.mammal)$Phylum,global_phy_cols$Phylum)]
### Assign Edge Colours Based on Positive Or Negativ
mammal_edgecols<-extract_edge_cols(mammals,mammal_speic,ig.mammal)
E(ig.mammal)$color=mammal_edgecols
#E(ig.mammal)$weight<-3
#Plot
plot(ig.mammal,vertex.size=9,vertex.label=NA)
plot_network(ig.mammal)
### Birds
birds <- subset_samples(microbial_merge,Taxa=="Bird") %>% filter_taxa(function(x) sum(x > 3) > (0.05*length(x)), TRUE)
bird_speic<- birds %>% spiec.easi(method='mb', lambda.min.ratio=1e-2, nlambda=20,pulsar.params=list(rep.num=50))
## Applying data transformations...
## Selecting model with pulsar using stars...
## Fitting final estimate with mb...
## done
ig.bird <- adj2igraph(getRefit(bird_speic),vertex.attr=list(name=taxa_names(birds),Kingdom=tax_table(birds)[,1],Phylum=tax_table(birds)[,2]))
#Assign Shapes to Kingdoms + Colours to Phylum
V(ig.bird)$shape<-ifelse(V(ig.bird)$Kingdom=="Bacteria","circle","square")
V(ig.bird)$color<-global_phy_cols$color[match(V(ig.bird)$Phylum,global_phy_cols$Phylum)]
### Assign Edge Colours Based on Positive Or Negativ
bird_edgecols<-extract_edge_cols(birds,bird_speic,ig.bird)
E(ig.bird)$color=bird_edgecols
#Plot
plot(ig.bird,vertex.size=9,vertex.label=NA,circular=T)
### Amphibians
amphibians <- subset_samples(microbial_merge,Taxa=="Amphibian") %>% filter_taxa(function(x) sum(x > 3) > (0.05*length(x)), TRUE)
amphib_speic<- amphibians %>% spiec.easi(method='mb', lambda.min.ratio=1e-2, nlambda=20,pulsar.params=list(rep.num=50))
## Applying data transformations...
## Selecting model with pulsar using stars...
## Fitting final estimate with mb...
## done
ig.amphib <- adj2igraph(getRefit(amphib_speic),vertex.attr=list(name=taxa_names(amphibians),Kingdom=tax_table(amphibians)[,1],Phylum=tax_table(amphibians)[,2]))
#Assign Shapes to Kingdoms + Colours to Phylum
V(ig.amphib)$shape<-ifelse(V(ig.amphib)$Kingdom=="Bacteria","circle","square")
V(ig.amphib)$color<-global_phy_cols$color[match(V(ig.amphib)$Phylum,global_phy_cols$Phylum)]
### Assign Edge Colours Based on Positive Or Negativ
amphibian_edgecols<-extract_edge_cols(amphibians,amphib_speic,ig.amphib)
E(ig.amphib)$color=amphibian_edgecols
#Plot
plot(ig.amphib,vertex.size=9,vertex.label=NA,main="Amphibians",circular=T)
#legend("topleft",legend=phylumcols_amphib[,1],fill=phylumcols_amphib[,2])
#### Fish
fish <- subset_samples(microbial_merge,Taxa=="Fish") %>% filter_taxa(function(x) sum(x > 3) > (0.05*length(x)), TRUE)
fish_speic<- fish %>% spiec.easi(method='mb', lambda.min.ratio=1e-2, nlambda=20,pulsar.params=list(rep.num=50))
## Applying data transformations...
## Selecting model with pulsar using stars...
## Fitting final estimate with mb...
## done
ig.fish <- adj2igraph(getRefit(fish_speic),vertex.attr=list(name=taxa_names(fish),Kingdom=tax_table(fish)[,1],Phylum=tax_table(fish)[,2]))
#Assign Shapes to Kingdoms + Colours to Phylum
V(ig.fish)$shape<-ifelse(V(ig.fish)$Kingdom=="Bacteria","circle","square")
V(ig.fish)$color<-global_phy_cols$color[match(V(ig.fish)$Phylum,global_phy_cols$Phylum)]
### Assign Edge Colours Based on Positive Or Negativ
fish_edgecols<-extract_edge_cols(fish,fish_speic,ig.fish)
E(ig.fish)$color=fish_edgecols
#Plot
plot(ig.fish,vertex.size=9,vertex.label=NA,main="Fish",circular=T)
#legend("topleft",legend=phylumcols_amphib[,1],fill=phylumcols_amphib[,2])
#### Insecta
sample_data(microbial_merge)$Class<-metadata$Class[match(sample_data(microbial_merge)$Common_name,metadata$Common_name)]
inverts <- subset_samples(microbial_merge,Class=="Insecta") %>% filter_taxa(function(x) sum(x > 3) > (0.05*length(x)), TRUE)
invert_speic<- inverts %>% spiec.easi(method='mb', lambda.min.ratio=1e-2, nlambda=20,pulsar.params=list(rep.num=50))
## Applying data transformations...
## Selecting model with pulsar using stars...
## Fitting final estimate with mb...
## done
ig.invert <- adj2igraph(getRefit(invert_speic),vertex.attr=list(name=taxa_names(inverts),Kingdom=tax_table(inverts)[,1],Phylum=tax_table(inverts)[,2]))
#Assign Shapes to Kingdoms + Colours to Phylum
V(ig.invert)$shape<-ifelse(V(ig.invert)$Kingdom=="Bacteria","circle","square")
V(ig.invert)$color<-global_phy_cols$color[match(V(ig.invert)$Phylum,global_phy_cols$Phylum)]
### Assign Edge Colours Based on Positive Or Negativ
invert_edgecols<-extract_edge_cols(inverts,invert_speic,ig.invert)
E(ig.invert)$color=invert_edgecols
#Plot
plot(ig.invert,vertex.size=9,vertex.label=NA,main="Insecta")
plot(ig.invert,vertex.size=9,vertex.label=NA,main="Invertebrate",circular=TRUE)
# jpeg(filename="network_plots.jpeg", width=300, height=320, units="mm", res=600)
final.fig <- par(mfcol=c(3, 2),mar=c(1,1,1,1),oma=c(1,1,1,1))
plot(ig.mammal,vertex.size=12,vertex.label=NA,main="Mammalia",circular=T,label.cex=10)
plot(ig.bird,vertex.size=12,vertex.label=NA,circular=T,main="Aves",circular=T)
plot(ig.amphib,vertex.size=12,vertex.label=NA,main="Amphibia",circular=T)
plot(ig.fish,vertex.size=12,vertex.label=NA,main="Actinopterygii",circular=T)
plot(ig.invert,vertex.size=12,vertex.label=NA,main="Insecta",circular=TRUE)
plot.new()
legend("topleft",legend=global_phy_cols[,1],fill=global_phy_cols[,2],border=NULL,cex=1.8,horiz=FALSE,ncol=2)
legend("bottomleft",legend=c("Bacteria","Fungi"),pch=c(21,22),cex=2.5,pt.cex=3)
# par(final.fig)
# dev.off()
#Assemble Data
taxa_edge_data<-data.frame(cols=c(mammal_edgecols,bird_edgecols,amphibian_edgecols,fish_edgecols,invert_edgecols),taxa=rep(c("Mammal","Bird","Amphibian","Fish","Invert"),times=c(length(mammal_edgecols),length(bird_edgecols),length(amphibian_edgecols),length(fish_edgecols),length(invert_edgecols))))
taxa_edge_data<- taxa_edge_data %>% mutate(sign=ifelse(cols=="forestgreen","Positive","Negative"))
taxa_edge_data %>% group_by(taxa) %>% summarise(npos=length(which(sign=="Positive")),nneg=length(which(sign=="Negative")))
## # A tibble: 5 x 3
## taxa npos nneg
## <chr> <int> <int>
## 1 Amphibian 68 10
## 2 Bird 241 6
## 3 Fish 64 1
## 4 Invert 121 0
## 5 Mammal 686 9
library(igraph)
#FUnction to Pull Out Edges and Filter to Bacteria-Fungi Interactions
edge_extract<-function(igraphobject){
#Pull Out Edge and Vertex Data
edges_data_frame <- get.data.frame(igraphobject, what = "edges")
v_data_frame <- get.data.frame(igraphobject, what = "vertices")
#Annotate With Metadata
edges_data_frame$king1<-v_data_frame$Kingdom[match(edges_data_frame$from,v_data_frame$name)]
edges_data_frame$king2<-v_data_frame$Kingdom[match(edges_data_frame$to,v_data_frame$name)]
edges_data_frame$kingpaste<-with(edges_data_frame,paste(king1,king2))
edges_data_frame$phylum1<-v_data_frame$Phylum[match(edges_data_frame$from,v_data_frame$name)]
edges_data_frame$phylum2<-v_data_frame$Phylum[match(edges_data_frame$to,v_data_frame$name)]
#Filter to Postive Co-Occurence
edges_bf<-edges_data_frame %>% filter(kingpaste=="Bacteria k__Fungi" & color=="forestgreen") %>% filter(phylum2!="Kingdom_Fungi")
return(edges_bf)
}
#### FUnction to Tabulate and Provide a Base Heatmap
heatplot<-function(edge_extract_output){
#Get Data
d1<-with(edge_extract_output,table(phylum1,phylum2)) %>% as.data.frame() %>% mutate(Prop=Freq/sum(Freq))
p1<- d1 %>% ggplot(aes(phylum1,phylum2,fill=Prop)) + geom_tile() + labs(x="Bacterial Phylum",y="Fungal Phylum") + plotopts + theme(axis.text.x = element_text(angle=45,hjust = 1))
return(p1)
}
#### Function to Permute Frequency of Co-occurrences
occur_rand<-function(edge_extract_output,nsims){
#Work Out Which is The Most Common
d1<-with(edge_extract_output,table(phylum1,phylum2)) %>% as.data.frame() %>% arrange(desc(Freq))
#Permutations
#Empty Vecgtor for Frequency
permute_vec<-numeric(nsims)
for(k in 1:nsims){
#Shuffle Phylum 2
edge_extract_output$phylum2shuffle<-sample(edge_extract_output$phylum2)
d_shuffle<-with(edge_extract_output,table(phylum1,phylum2shuffle)) %>% as.data.frame()
permute_vec[k]<-d_shuffle$Freq[which(d_shuffle$phylum1==d1[1,1] & d_shuffle$phylum2shuffle==d1[1,2])]
} #loop end
print(paste("The most common co-ccurence was",d1[1,1],"&",d1[1,2]))
twosidep<-mean(permute_vec<= min(d1[1,3],-d1[1,3])) + mean(permute_vec>= max(d1[1,3],-d1[1,3]))
print(paste(" 2 tailed p value = ",twosidep))
} #function end
#### MAMMALS
mammal_bf<-edge_extract(ig.mammal)
heatplot(mammal_bf) + scale_fill_gradient()
occur_rand(mammal_bf,10000)
## [1] "The most common co-ccurence was Firmicutes & Ascomycota"
## [1] " 2 tailed p value = 0.7595"
#### BIRDS
bird_bf<-edge_extract(ig.bird)
heatplot(bird_bf) + scale_fill_gradient()
occur_rand(bird_bf,10000)
## [1] "The most common co-ccurence was Actinobacteria & Ascomycota"
## [1] " 2 tailed p value = 0.0014"
#### INVERTS
invert_bf<-edge_extract(ig.invert)
heatplot(invert_bf) + scale_fill_gradient()
occur_rand(invert_bf,10000)
## [1] "The most common co-ccurence was Proteobacteria & Ascomycota"
## [1] " 2 tailed p value = 0.6201"
#### FISH
fish_bf<-edge_extract(ig.fish)
heatplot(fish_bf) + scale_fill_gradient()
occur_rand(fish_bf,10000)
## [1] "The most common co-ccurence was Proteobacteria & Basidiomycota"
## [1] " 2 tailed p value = 1"
### AMPHIBIA
amphib_bf<-edge_extract(ig.amphib)
heatplot(amphib_bf) + scale_fill_gradient()
occur_rand(amphib_bf,10000)
## [1] "The most common co-ccurence was Proteobacteria & Ascomycota"
## [1] " 2 tailed p value = 1"
## Adjust P
p.adjust(c(0.765,0.002,0.6205,1,1))
## [1] 1.00 0.01 1.00 1.00 1.00
mammal_heat<-heatplot(mammal_bf) + scale_fill_gradient() + labs(title="Mammalia",fill="Proportion") +
theme(plot.title = element_text(hjust = 0.5,face="italic",size=15),legend.title = element_text(size=12))
bird_heat<-heatplot(bird_bf) + scale_fill_gradient()+ labs(title="Aves",fill="Proportion") +
theme(plot.title = element_text(hjust = 0.5,face="italic",size=15),legend.title = element_text(size=12))
invert_heat<-heatplot(invert_bf) + scale_fill_gradient()+ labs(title="Insecta",fill="Proportion") +
theme(plot.title = element_text(hjust = 0.5,face="italic",size=15),legend.title = element_text(size=12))
fish_heat<-heatplot(fish_bf) + scale_fill_gradient()+ labs(title="Actinopterygii",fill="Proportion") +
theme(plot.title = element_text(hjust = 0.5,face="italic",size=15),legend.title = element_text(size=12))
amphibian_heat<-heatplot(amphib_bf) + scale_fill_gradient()+ labs(title="Amphibia",fill="Proportion") +
theme(plot.title = element_text(hjust = 0.5,face="italic",size=15),legend.title = element_text(size=12))
heat_grid<-plot_grid(mammal_heat,bird_heat,invert_heat,fish_heat,amphibian_heat,nrow=3,ncol=2)
heat_grid
ggsave2('Heatmap Grid.pdf',heat_grid,width=10,height=10,units="in")
#Read in Co-oCcur Data
posneg<-read.csv('posneg.csv',header=T)
head(posneg)
## Common_name Positive Negative
## 1 African palm weevil larvae 14 9
## 2 Bank vole 43 25
## 3 Barred hamlet 21 24
## 4 Blackcap 188 231
## 5 Capuchin monkey 78 64
## 6 Carrion crow 31 35
posneg$Common_name<-as.character(posneg$Common_name)
#posneg$Class<-metadata$Class[match(posneg$Common_name,metadata$Common_name)]
metadata$Common_name<-as.character(metadata$Common_name)
#Combine With Metadata
head(metadata)
## Sample2 Species Common_name Taxa Sex Situ
## 1 Alytes-1 Alytes obstetricans Common midwife toad Amphibian Unknown Captive
## 2 Alytes-13 Alytes obstetricans Common midwife toad Amphibian Unknown Captive
## 3 Alytes-14 Alytes obstetricans Common midwife toad Amphibian Unknown Captive
## 4 Alytes-15 Alytes obstetricans Common midwife toad Amphibian Unknown Captive
## 5 Alytes-16 Alytes obstetricans Common midwife toad Amphibian Unknown Captive
## 6 Alytes-2 Alytes obstetricans Common midwife toad Amphibian Unknown Captive
## Location Country CollectionYear TissueStorage SampleType
## 1 London Zoo, UK Captive 2015 Frozen at -20C Skin swab
## 2 London Zoo, UK Captive 2015 Frozen at -20C Skin swab
## 3 London Zoo, UK Captive 2015 Frozen at -20C Skin swab
## 4 London Zoo, UK Captive 2015 Frozen at -20C Skin swab
## 5 London Zoo, UK Captive 2015 Frozen at -20C Skin swab
## 6 London Zoo, UK Captive 2015 Frozen at -20C Skin swab
## ExtractionKit Environment Class Order Rough_Class Diet
## 1 Qiagen DNEasy kit Freshwater Amphibia Anura Amphibia Carnivore
## 2 Qiagen DNEasy kit Freshwater Amphibia Anura Amphibia Carnivore
## 3 Qiagen DNEasy kit Freshwater Amphibia Anura Amphibia Carnivore
## 4 Qiagen DNEasy kit Freshwater Amphibia Anura Amphibia Carnivore
## 5 Qiagen DNEasy kit Freshwater Amphibia Anura Amphibia Carnivore
## 6 Qiagen DNEasy kit Freshwater Amphibia Anura Amphibia Carnivore
posnegj<- metadata %>% select("Common_name","Taxa","Class","Order") %>% unique %>% inner_join(x=posneg,"Common_name")
posnegj<- mutate(posnegj,n_interactions=Positive + Negative)
#### MODELS - Binomial
posnegj_filter <- posnegj %>% filter(Class %in% c("Actinopterygii","Aves","Mammalia","Insecta","Amphibia"))
degree_prop1<- brm(Positive|trials(n_interactions) ~ Class,family=binomial(),data=posnegj_filter)
## Compiling Stan program...
## Start sampling
summary(degree_prop1)
## Family: binomial
## Links: mu = logit
## Formula: Positive | trials(n_interactions) ~ Class
## Data: posnegj_filter (Number of observations: 39)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 0.46 0.17 0.11 0.80 1.01 774 1228
## ClassAmphibia -0.20 0.20 -0.59 0.21 1.01 950 1789
## ClassAves -0.13 0.18 -0.48 0.22 1.01 784 1189
## ClassInsecta -0.60 0.23 -1.06 -0.16 1.00 1029 1843
## ClassMammalia 0.05 0.18 -0.30 0.39 1.01 813 1551
##
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
pp_check(degree_prop1)
## Using 10 posterior samples for ppc type 'dens_overlay' by default.
conditional_effects(degree_prop1)
## Setting all 'trials' variables to 1 by default if not specified otherwise.
bayes_R2(degree_prop1)
## Estimate Est.Error Q2.5 Q97.5
## R2 0.9314302 0.001351038 0.9284625 0.9336376
## Taxa Names
taxlabels<-c("Actinopterygii","Amphibia","Aves","Insecta","Mammalia")
########## Binomial Model
propmod_plotdata<-fixef(degree_prop1,summary=F)
propmod_plotdata2<-data.frame(plogis(cbind(propmod_plotdata[,1],propmod_plotdata[,1]+propmod_plotdata[,2:5])))
colnames(propmod_plotdata2)<-taxlabels
propmod_plotdata3<- propmod_plotdata2 %>% pivot_longer(colnames(propmod_plotdata2),names_to="Class",values_to="value")
propmod_plotdata3$facetval<-"proportion"
###Taxonomic Diffs
# bird_mamm_prop<-propmod_plotdata2$Bird - propmod_plotdata2$Mammal
# mean(bird_mamm_prop); quantile(bird_mamm_prop,c(0.025,0.975))
# mammal_amphib_prop<-propmod_plotdata2$Mammal - propmod_plotdata2$Amphibian
# mean(mammal_amphib_prop); quantile(mammal_amphib_prop,c(0.025,0.975))
#
#
######### Plot - Proportion Positive Edges
network_statplot2<- ggplot(propmod_plotdata3,aes(x=value,y=Class)) + geom_vline(xintercept = 0.5,lty="dashed") + theme_bw() + plotopts + stat_halfeye(shape=21,point_fill="white",point_size=5,colour="black",slab_alpha=0.5,aes(fill=Class)) + guides(fill=F)+
labs(x="Proportion of Positive Edge Correlations")
#ggsave2('Netowrk Edge Proportion Graph.pdf',network_statplot2,width=15,height=15,units="cm")
# ggsave2('Netowrk Edge Proportion Graph.jpeg',network_statplot2,width=15,height=15,units="cm")
#Estimate Diversity - Fungi
network_rich_fungi<- myco %>% subset_samples(Common_name %in% posnegj$Common_name) %>% rarefy_even_depth(sample.size=500) %>% estimate_richness(measures = c("Shannon"))
## You set `rngseed` to FALSE. Make sure you've set & recorded
## the random seed of your session for reproducibility.
## See `?set.seed`
## ...
## 51 samples removedbecause they contained fewer reads than `sample.size`.
## Up to first five removed samples are:
## Alytes-1Alytes-13Alytes-15Alytes-16Alytes-2
## ...
## 67OTUs were removed because they are no longer
## present in any sample after random subsampling
## ...
network_rich_fungi$Common_name<- metadata$Common_name[match(rownames(network_rich_fungi),metadata$Sample2)]
#Estimate Diversity - Bactera
network_rich_bacteria<- micro %>% subset_samples(Common_name %in% posnegj$Common_name) %>% rarefy_even_depth(sample.size=500) %>% estimate_richness(measures = c("Shannon"))
## You set `rngseed` to FALSE. Make sure you've set & recorded
## the random seed of your session for reproducibility.
## See `?set.seed`
##
## ...
## 26 samples removedbecause they contained fewer reads than `sample.size`.
## Up to first five removed samples are:
## Alytes-15Alytes-16Bird25Bird29Bird62
## ...
## 17OTUs were removed because they are no longer
## present in any sample after random subsampling
## ...
network_rich_bacteria$Common_name<- metadata$Common_name[match(rownames(network_rich_bacteria),metadata$Sample2)]
#Mean Values
network_rich_fungi_mean<- network_rich_fungi %>% group_by(Common_name) %>% summarise(mean_shannon=mean(Shannon))
network_rich_bacteria_mean<- network_rich_bacteria %>% group_by(Common_name) %>% summarise(mean_shannon=mean(Shannon))
#Copy Across
posnegj$mean_fungal_shannon<-network_rich_fungi_mean$mean_shannon[match(posnegj$Common_name,network_rich_fungi_mean$Common_name)]
posnegj$mean_bacteria_shannon<-network_rich_bacteria_mean$mean_shannon[match(posnegj$Common_name,network_rich_bacteria_mean$Common_name)]
head(posnegj)
## Common_name Positive Negative Taxa Class
## 1 African palm weevil larvae 14 9 Invertebrate Insecta
## 2 Bank vole 43 25 Mammal Mammalia
## 3 Barred hamlet 21 24 Fish Actinopterygii
## 4 Blackcap 188 231 Bird Aves
## 5 Capuchin monkey 78 64 Mammal Mammalia
## 6 Carrion crow 31 35 Bird Aves
## Order n_interactions mean_fungal_shannon mean_bacteria_shannon
## 1 Coleoptera 23 0.3933196 2.3824486
## 2 Rodentia 68 0.6383061 2.4488759
## 3 Perciformes 45 0.7792424 1.8451177
## 4 Passeriformes 419 0.9932420 1.1894112
## 5 Primates 142 0.9083341 1.4398882
## 6 Passeriformes 66 0.6840271 0.6922338
# Calculate Mean Proportion
posnegj$proportion<- with(posnegj,Positive / n_interactions)
#Filter to The Big 5
posnegj_filter <- posnegj %>% filter(Class %in% c("Actinopterygii","Aves","Mammalia","Insecta","Amphibia"))
######### PLOTS - SPECIES LEVEL
proportion_shannon_bacteria_plot1<- ggplot(posnegj_filter,aes(x=mean_bacteria_shannon,y=proportion)) + geom_point(aes(color=Class))
with(posnegj_filter,cor.test(mean_bacteria_shannon,proportion))
##
## Pearson's product-moment correlation
##
## data: mean_bacteria_shannon and proportion
## t = -0.35308, df = 37, p-value = 0.726
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.3667596 0.2623658
## sample estimates:
## cor
## -0.05794821
proportion_shannon_fungi_plot1<- ggplot(posnegj_filter,aes(x=mean_fungal_shannon,y=proportion)) + geom_point()
with(posnegj_filter,cor.test(mean_fungal_shannon,proportion))
##
## Pearson's product-moment correlation
##
## data: mean_fungal_shannon and proportion
## t = 0.31016, df = 37, p-value = 0.7582
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.2689142 0.3606460
## sample estimates:
## cor
## 0.05092353
######### PLOTS - CLASS LEVEL
posneg_class<- posnegj_filter %>% group_by(Class) %>% summarise(mean_prop=mean(proportion),mean_bac_shannon=mean(mean_bacteria_shannon),mean_fungal_shannon=mean(mean_fungal_shannon))
proportion_shannon_bacteria_classplot1<- ggplot(posneg_class,aes(x=mean_bac_shannon,y=mean_prop)) + geom_smooth(method="lm") + geom_point(shape=21,size=4,aes(fill=Class)) + labs(x="Mean Bacterial Shannon Diversity",y="Mean Proportion Positive Edges") + plotopts
proportion_shannon_bacteria_classplot1
## `geom_smooth()` using formula 'y ~ x'
with(posneg_class,cor.test(mean_bac_shannon,mean_prop))
##
## Pearson's product-moment correlation
##
## data: mean_bac_shannon and mean_prop
## t = 0.2051, df = 3, p-value = 0.8506
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.8531904 0.9058763
## sample estimates:
## cor
## 0.1175935
proportion_shannon_fungi_classplot1<- ggplot(posneg_class,aes(x=mean_fungal_shannon,y=mean_prop)) + geom_smooth(method="lm") + geom_point(shape=21,size=4,aes(fill=Class))+ labs(x="Mean Fungal Shannon Diversity",y="Mean Proportion Positive Edges") + plotopts
proportion_shannon_fungi_classplot1
## `geom_smooth()` using formula 'y ~ x'
with(posneg_class,cor.test(mean_fungal_shannon,mean_prop))
##
## Pearson's product-moment correlation
##
## data: mean_fungal_shannon and mean_prop
## t = 2.1157, df = 3, p-value = 0.1247
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.3418874 0.9841717
## sample estimates:
## cor
## 0.7737782
Here we randomly reassign node labels across kingdoms (bacteria & fungi) to test whether fungi are critical betweenness components
######### RANDOMISATIONS
#Number of Shuffles
nrand<-1000
###########MAMMALS
###Betweenness of Fungi
mammal_between<-betweenness(ig.mammal)
mammal_fungal_permute<- replicate(10000,mean(mammal_between[grep("fungi",sample(names(mammal_between)))]))
mammal_b_true<- mean(mammal_between[grep("fungi",names(mammal_between))])
mean(mammal_fungal_permute>=mammal_b_true)*2
## [1] 0.6292
#Random Removal of Bacteria or Fungi
mammal_nfungi<-length(which(tax_table(mammals)[,1]=="k__Fungi"))
sample(seq(nrow(tax_table(mammals))),mammal_nfungi)
## [1] 100 128 215 90 38 74 217 11 29 34 117 133 25 226 184 172 141 77 224
## [20] 20 4 126 156 169 137 183 44 101 58 132 10 160
mammal_no_fungi<- delete_vertices(ig.mammal,which(tax_table(mammals)[,1]=="k__Fungi")) %>% cluster_fast_greedy() %>% modularity()
qr_mammal_fungidelete_perm <- replicate(10000, delete_vertices(ig.mammal,sample(seq(nrow(tax_table(mammals))),mammal_nfungi)),simplify = F) %>% lapply(cluster_fast_greedy) %>%
sapply(modularity)
mean(mammal_no_fungi>=qr_mammal_fungidelete_perm) #0.0441
## [1] 0.0464
##Betweenness of Fungi
mammal_bf_bacbetween<- delete_vertices(ig.mammal,which(tax_table(mammals)[,1]=="k__Fungi")) %>% betweenness() %>% mean()
qr_mammal_fungidelete_between_perm <- replicate(1000, delete_vertices(ig.mammal,sample(seq(nrow(tax_table(mammals))),mammal_nfungi)),simplify = F) %>% sapply(betweenness) %>% apply(2,mean)
hist(qr_mammal_fungidelete_between_perm)
mean(mammal_bf_bacbetween<=qr_mammal_fungidelete_between_perm)
## [1] 0.144
###########BIRDS
#Betweenness of Fungi
bird_between<-betweenness(ig.bird)
bird_fungal_permute<- replicate(10000,mean(bird_between[grep("fungi",sample(names(bird_between)))]))
bird_b_true<- mean(bird_between[grep("fungi",names(bird_between))])
mean(bird_fungal_permute<=bird_b_true)*2 #0.04
## [1] 0.047
#Hist
hist(bird_fungal_permute)
#Assemble Data
between_permute_data<-data.frame(randval=c(mammal_fungal_permute,bird_fungal_permute),Class=rep(c("Mammalia","Aves"),each=nrand))
# #Plots
mammal_permplot1<- between_permute_data %>% filter(Class=="Mammalia") %>% ggplot(.,aes(randval)) + geom_density(fill="lightgray")
mammal_permplot2<- mammal_permplot1<- mammal_permplot1 + theme_bw() + plotopts + labs(x="Mammal Randomised Network Betweenness",y="Density")
mammal_permplot3<- mammal_permplot2 + geom_vline(xintercept = mammal_b_true,cex=1.5,col="red")
bird_permplot1<- between_permute_data %>% filter(Class=="Aves") %>% ggplot(.,aes(randval)) + geom_density(fill="lightgray")
bird_permplot2<- bird_permplot1 + theme_bw() + plotopts + labs(x="Bird Randomised Network Betweenness",y="Density")
bird_permplot3<- bird_permplot2 + geom_vline(xintercept = bird_b_true,cex=1.5,col="red")
#Combine
dummydf<-data.frame(Class=c("Mammalia","Aves"),realval=c(mammal_b_true,bird_b_true)) #real value
textdf<-data.frame(label= c("p = 0.044*","p = 0.63"),Class=c("Aves","Mammalia"),x=750,y=0.003)
#textdf2<-data.frame(label= c("2-tailed p = 0.044","2-tailed p = 0.63"),Class=c("Aves","Mammalia"),x=850,y=0.003)
perm_plot_all<-plot_grid(mammal_permplot3,bird_permplot3,labels=c("Mammalia","Aves"),ncol=2)
permplot1<-ggplot(between_permute_data,aes(randval)) + geom_density(aes(fill=Class),alpha=0.5) + facet_wrap(.~Class)
permplot2<- permplot1 + theme_bw() + plotopts + guides(fill=F) + labs(x="Randomised Betweenness Value",y="Density")
permplot3<- permplot2 + geom_vline(data=dummydf,aes(xintercept=realval),cex=1.5,color="red") + geom_text(data = textdf,
mapping = aes(x = x,
y = y,
label = label),size=8)
### Function That Strips Out Fungal Nodes for Betweenness Calculations
fungus_subset<-function(matrix_in,physeq_obj){
return(matrix_in[tax_table(physeq_obj)[,1]=="k__Fungi",])
}
### Function That Calculates True Fungal Betweenness
extract_fungal_between<-function(igraph_obj,physeq_obj){
return(betweenness(igraph_obj)[tax_table(physeq_obj)[,1]=="k__Fungi"])
}
############# Modularity and Group Calcs
modularity_data<-data.frame(class=c("Mammalia","Aves","Insecta","Actinopterygii","Amphibia"))
#Cluster Each IGraph PObject
mammal_cluster<-ig.mammal %>% cluster_fast_greedy()
bird_cluster<-ig.bird %>% cluster_fast_greedy()
insect_cluster<-ig.invert %>% cluster_fast_greedy()
fish_cluster<-ig.fish %>% cluster_fast_greedy()
amphib_cluster<-ig.amphib %>% cluster_fast_greedy()
#Strip Out Data
modularity_temp<-c(modularity(mammal_cluster),modularity(bird_cluster),modularity(insect_cluster),modularity(fish_cluster),modularity(amphib_cluster))
groups_temp<-c(max(mammal_cluster$membership),max(bird_cluster$membership),max(insect_cluster$membership),max(fish_cluster$membership),max(amphib_cluster$membership))
components_temp<-c(components(ig.mammal)$no,components(ig.bird)$no,components(ig.invert)$no,components(ig.fish)$no,components(ig.amphib)$no)
#Add TO Data Frame
modularity_data$modularity<-round(modularity_temp,3)
modularity_data$groups<-groups_temp
modularity_data$components<-components_temp
#Stats
with(modularity_data,cor(modularity,groups))
## [1] 0.3947702
with(modularity_data,cor(modularity,components))
## [1] 0.6737694
#Plot
mod_groups_plot<-ggplot(modularity_data,aes(x=groups,y=modularity)) + geom_smooth(method="lm") + geom_point(size=5,shape=21,fill="white") + theme_bw() + plotopts + labs(x="Number of Groups",y="Modularity")
mod_components_plot<-ggplot(modularity_data,aes(x=components,y=modularity)) + geom_smooth(method="lm") + geom_point(size=5,shape=21,fill="white") + theme_bw() + plotopts + labs(x="Number of Components",y="Modularity")
mod_combined_plot<-plot_grid(mod_groups_plot,mod_components_plot,labels=c("A","B"),label_size=30)
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
ggsave2('Modularity Combined Plot.pdf',mod_combined_plot,width=35,height=15,units="cm")
mod_combined_plot
###########Modularity Distributions
mammal_modularity_dist<-replicate(10000, delete_vertices(ig.mammal,sample(seq(nrow(tax_table(mammals))),round(length(V(ig.mammal))*0.1))),simplify = F) %>% lapply(cluster_fast_greedy) %>%
sapply(modularity)
bird_modularity_dist<-replicate(10000, delete_vertices(ig.bird,sample(seq(nrow(tax_table(birds))),round(length(V(ig.bird))*0.1))),simplify = F) %>% lapply(cluster_fast_greedy) %>%
sapply(modularity)
invert_modularity_dist<-replicate(10000, delete_vertices(ig.invert,sample(seq(nrow(tax_table(inverts))),round(length(V(ig.invert))*0.1))),simplify = F) %>% lapply(cluster_fast_greedy) %>%
sapply(modularity)
fish_modularity_dist<-replicate(10000, delete_vertices(ig.fish,sample(seq(nrow(tax_table(fish))),round(length(V(ig.fish))*0.1))),simplify = F) %>% lapply(cluster_fast_greedy) %>%
sapply(modularity)
amphib_modularity_dist<-replicate(10000, delete_vertices(ig.amphib,sample(seq(nrow(tax_table(amphibians))),round(length(V(ig.amphib))*0.1))),simplify = F) %>% lapply(cluster_fast_greedy) %>%
sapply(modularity)
######## Assemble Modularity Data
modularity_plotdata<-data.frame(modularity=c(mammal_modularity_dist,bird_modularity_dist,invert_modularity_dist,fish_modularity_dist,amphib_modularity_dist),class=rep(c("Mammalia","Aves","Insecta","Actinopterygii","Amphibia"),each=10000))
######## Plot
ggplot(modularity_plotdata,aes(x=modularity,y=class)) + stat_halfeye(shape=21,point_fill="white",point_size=5,colour="black",slab_alpha=0.5,aes(fill=class))
######### Plot - Modularity
network_statplot1<- ggplot(modularity_plotdata,aes(x=modularity,y=class)) + theme_bw() + plotopts + stat_halfeye(shape=21,point_fill="white",point_size=5,colour="black",slab_alpha=0.5,aes(fill=class)) + guides(fill=F)+ labs(x="Modularity",y="Class")
######### Plot - Proportion Positive Edges
network_statplot2<- ggplot(propmod_plotdata3,aes(x=value,y=Class)) + geom_vline(xintercept = 0.5,lty="dashed") + theme_bw() + plotopts + stat_halfeye(shape=21,point_fill="white",point_size=5,colour="black",slab_alpha=0.5,aes(fill=Class)) + guides(fill=F)+
labs(x="Proportion of Positive Edge Correlations")
#ggsave2('Netowrk Edge Proportion Graph.pdf',network_statplot2,width=15,height=15,units="cm")
# ggsave2('Netowrk Edge Proportion Graph.jpeg',network_statplot2,width=15,height=15,units="cm")
#Extract Colour Data fro Classes from Proportion Neg Plot
g <- ggplot_build(network_statplot2)
class_fill_cols<-unique(g$data[[2]]["fill"])[,1]
#Rejig Permplot with custom fill
classcols<-c("Mammalia"=class_fill_cols[5],"Aves"=class_fill_cols[3])
permplot4<- permplot3 + scale_fill_manual(values=classcols)
## Part 1 - Subplot of Modularity, Betweeness Randomisations and Proportion Positive Edges
network_plot_sub1<-plot_grid(network_statplot2,network_statplot1,permplot4,labels=c("A","B","C"),label_size = 30,ncol=1,nrow=3)
#Part 2 - Add Heatmap
network_plot_sub2<- plot_grid(network_plot_sub1,heat_grid,ncol=2,labels=c("","D"),label_size=30,rel_widths = c(1.4,2))
network_plot_sub2
#Save
ggsave2('Combined Network Statistics Graph v4 Aug21.pdf',network_plot_sub2,width=42,height=30,units="cm")
ggsave2('Combined Network Statistics Graph v4 Jul21.jpeg',network_plot_sub2,width=18,height=12,units="in")
ggsave2('Combined Network Statistics Graph v4 Jul21.tiff',network_plot_sub2,width=18,height=12,units="in")